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ABSTRACT 

The  unsteady,  viscous,  three-dimensional  flow  field  of  a  he¬ 
licopter  rotor  blade  encountering  a  passing  vortex  is  calculated 
by  solving  the  thin  layer  Navier-Stokes  equations  by  a  finite- 
difference  numerical  procedure.  A  prescribed  vortex  method  is 
adopted  to  preserve  the  structure  q{  the  interacting  concentrated 
vortex.  The  test  cases  considered  correspond  to  the  experimen¬ 
tal  model  rotor  test  condition^  of  Caradonnaet  al.’’^ 

Uf,  .  's,  '' 

NOMENCLATURE 


Ugo 

= 

characteristic  velocity  scale,  speed  of  sound 

C 

= 

characteristic  length  scale,  chord  of  the  blade 

C, 

= 

pressure  coefficient 

e 

= 

total  energy  per  unit  volume 

F,G,H 

flux  ye^tors 

free-stream  Mach  number 

= 

tip  Mach  number  of  the  rotor  blade 

P 

= 

static  pressure 

Pr 

= 

Prandtl  number 

Q 

= 

flowfield  vector 

'■fl 

= 

rotor  reference  station 

n 

= 

rotor  radius 

mi) 

= 

rotational  matrix,  see  Eq.  (7) 

Re 

= 

Reynolds  number 

s 

= 

viscous  flux  vector 

= 

free  stream  velocity 

u,»,u> 

= 

velocity  components 

v,v,w 

= 

contravariant  velocity  components 

= 

inertial  coordinates 

= 

blade  attached  coordinates 

ratio  of  specific  heats 
dimensionless  strength  of  tip  vortex 
pitch  angle  of  the  rotor  blade 
viscosity  coefficient,  advance  ratio 
generalized  curvilinear  coordinates 
density 
azimuth  angle 
angular  velocity  of  the  rotor 


INTRODUCTION 

The  accurate  simulation  of  the  flowfield  of  a  helicopter  ro¬ 
tor  is  still  one  of  the  most  complex  and  challenging  problems 
of  applied  aerodynsonics.  This  is  true  inspite  of  the  availabil¬ 
ity  of  the  present  day  supercomputers  of  Cray-2  class  and  im¬ 
proved  numerical  algorithms.  The  flowfield  is  three-dimensional, 
unsteady  and  viscous,  with  pockets  of  transonic  flow  near  the 
blade  tips  on  the  advancing  side  and  regions  of  i .  namic  stall 
pockets  on  the  retreating  blades.  In  addition,'  the  blades  also 
shed  complex  vortical  wakes;  the  concentrated  tip  vortex  of  each 
blade  generally  passes  close  to  the  following  blades.  The  close 
encounter  of  these  (force  free)  concentrated  vortices  with  the  ro¬ 
tor  blades  is  often  the  cause  of  unsteady  load  fluctuations  and 
acoustic  noise.  Figure  1  shows,  through  a  schematic,  the  various 
complexities  associated  with  the  flowfield  of  a  typical  helicopter 
rotor  in  forward  flight. 

The  spiralling  vortex  sheet  emanating  from  each  of  the 
blades  of  the  rotor  has  a  profound  influence  on  the  performance 
of  the  helicopter.  It  not  only  alters  the  effective  pitch  angle  of  the 
blades  and  thus  the  airloads,  but  also  produces  highly  nonlinear 
interactions  of  the  vortex  with  the  rotor  flowfield.  This  might 
result  many  times  in  shock  induced  boundary  layer  separation 
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with  associated  loss  of  lift  and  increased  drag.  An  accurate 
simulation  of  the  rotor  flowfield,  therefore,  must  consider  the 


ir>icn 


Fig.  1  Schematic  of  complex  flowfield  of  a  helicopter  rotor 

induced  effects  of  the  vortex  wake  including  effects  of  blade- 
vortex  interactions  (BVI).  Numerical  simtdation  of  vortex  wakes 
is  being  attempted  only  recently  after  the  bigger  and  faster  com¬ 
puters  have  become  available  (see  for  example  Ref.  3)  and  have 
some  limited  success  so  far.  On  the  otherhand,  much  of  the 
progress  made  to  date  in  modelling  the  blade- vortex  interaction 
has  been  hampered  by  the  lack  of  the  development  of  theoreti¬ 
cal  and/or  numerical  techniques  to  preserve  the  structure  of  the 
concentrated  vortices  in  the  flowfield  without  significant  diffu¬ 
sion.  The  study  of  blade- vortex  interaction  has  been  the  subject 
of  numerous  research  papers  recently.^"*  Among  these  studies, 
some  of  the  methods  proposed’’~^  are  not  only  very  promising  in 
preserving  the  concentrated  vortices  but  have  actually  demon¬ 
strated  to  perform  satisfactorily,  atleast,  in  two-dimensional  in¬ 
teractions.  For  example,  Fig.  2  shows  a  comparison  of  experi¬ 
mental  data  of  Caradonna  et  al.*  with  the  unsteady,  viscous. 


two-dimensional  calculations  of  Srinivasan  et  al.*.  This  method 
used  a  prescribed  vortex  approach  for  preserving  the  convecting 
vortex.^  Although  for  the  subcritical  flow  condition  of  Fig.  2,  the 
comparison  of  calculations  with  data  is  excellent,  at  the  super¬ 
critical  flow  conditions  the  shock  wave  dominates  the  flowfield 
and  the  time  lag  effects  in  shock  wave  growth  and  subsequent 
decay  become  a  dominant  factor.  In  s^ldition,  three-dimensional 
effects  get  accentuated  in  the  presence  of  shock  waves.  Hence  a 
two-dimensional  representation  of  such  a  flowfield  is  inadequate 
as  observed  in  the  study  of  Ref.  8. 

The  purpose  of  this  paper  is  to  calculate  accurately  the  un¬ 
steady  interaction  flowfield  measured  in  a  laboratory  experiment 
by  Caradonna  et  al.’*^  on  a  model  helicopter  rotor  consisting  of 
both  the  parallel  and  skewed  blade-vortex  interactions.  In  this 
experiment,  the  interacting  vortex  was  generated  upstream  of 
the  model  rotor  by  means  of  a  wing. 

GOVERNING  EQUATIONS  AND  NUMERICAL  SCHEME 

The  unsteady,  viscous  flowfield  calculations  are  performed 
using  a  recently  developed  thin  layer  Navier-Stokes  code  specif¬ 
ically  written  for  helicopter  rotor  applications. The  govern¬ 
ing  partial  differential  equations  are  the  unsteady,  thin-layer 
Navier-Stokes  equations.  For  generality,  these  equations  are 
transformed  from  the  Cartesian  reference  frame  (  z,  y,  z,  f  )  to 
the  arbitrary  curvilinear  space  ((,  q,  C,  r)  while  retaining  strong 
conservation  law-form  to  capture  shock  waves.  The  transformed 
equations,  following  Ref.  7,  can  be  written  as 

dr{Q  -  $,)  +  ddF  -  F^)  +  -  G,)  +  -  H^) 

=  Re-^d(S 


CMOMAJiom  ixnMMfTTs  mKfita 

-  •  Lowm 

-  o// 


-  O.C.  /i  =  0.2. 


where 


■  p  ■ 

pU 

pu 

puU  -t-{,p 

pv 

pVU  +  (yP 

pw 

pxvU  +  (.p 

.  e  . 

.U[e  +  p)  -  i,p. 

pV  1 

1  1 

p\v 

puV  +  q,p 

1 

il 

puW  +  CrP 

G  =  J-' 

pvV  +  rjyp 

pvW  -l-<,p 

pwV  +  i/.p 

puilV  -f  C.p 

.V(e-(-p)-77,p. 

Lir(e  +  p)-(;.pJ 

Here  Q  is  the  flowfield  vector  we  are  solving  for  and  Q,  is  the 
buiulioii  of  the  Euler  equations  for  a  line  vortex  convecting  in  a 
uniform  free  stream.  F,  G,  11  and  G«,  are  the  appro¬ 
priate  flux  vectors  for  the  two  flows.  The  contravariant  velocity 
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components  U,  V  and  W  are  defined  as 
U  =  +  Cv"  + 

V  =  T]t  +  »7,U  +  TfyV  +  7]tW 

W  =  C<  +  <*«  +  C»»  +  C*® 


The  viscous  flux  vector  S,  written  here  in  the  limit  of  thin-layer 
approximation,  is  given  by 


are  related  to  the  blade  fixed  coordinates  Xb  =  {x,y,z,i) 
through  the  relation  given  by 

X(z,y,z)  =  R(t)Xb(i,y,z) 

(3) 

l  =  t  '  ' 


where  R(t)  is  the  rotational  matrix*^  given  by 


m  = 


cos  Clt 
sin  Of 
0 


— sin  ni  0 

cos  (It  0 
0  1 


(4) 


0 


S  =  J-‘ 


-I-  (/i/3)mj<* 
lxm\W(  -i-  {fi/3)mi^, 

./tmims  +  (/i/3)mj(<,u  +  C,u  +  C^). 


where 


mt  =  -I-  (yV(  + 

u*  -I-  +  KPr~^{y  -  l)~‘(a^)^ 

In  this  equation,  k  is  the  coefficient  of  thermal  conductivity,  Re 
is  the  Reynolds  number,  Pr  is  the  Prandtl  number,  a  is  the  speed 
of  sound  and  J  is  the  transformation  Jacobian,  whose  inverse  is 

+  X(y(Z„  +  x„y<2j  -  xty<z. 


The  primitive  variables  of  Eq.  (1),  viz.,  the  density  p,  the 
mass  fluxes  pu,pv,pw  and  the  energy  per  unit  volume  e,  are 
normalized  by  the  free  stream  reference  quantities.  The  refer¬ 
ence  length  scale  is  C  and  the  reference  time  scale  is  C/aga- 
The  details  on  the  m^Crics  of  the  transformation 

and  can  be  found  in  Ref.ill. 

The  velocity  components  u,  v,u;  and  the  pressure,  p,  are 
related  to  the  total  energy  per  unit  volume,  e,  through  the  equa¬ 
tion  of  state  for  a  perfect  gas  by 


p  =  (7 -l)(e-^(u’ +«*+«;’))  (2) 

where  7  is  the  ratio  of  specific  heats. 

In  the  above  equations  u,v,  and  w  are  the  Cartesian 
romponciils  of  llic  velocity  in  the  inertial  coordinate  system 
{r,y,z,t).  In  the  present  formulation  Eq.  (1)  is  solved  in  the  in¬ 
ertial  frame  of  reference.  The  inertial  coordinates  X  =  (i,  y,  z,  t) 


Here  fl  is  the  rotational  frequency  of  the  rotor  and  Rt  represents 
the  azimuth  sweep  of  the  rotor  blade. 

The  equation  set,  Eq.  (1),  is  solved  using  an  implicit, 
approximately-factored  numerical  scheme,  developed  by  Ying  et 
al.,**  that  uses  spatial  central-differencing  in  the  7  and  (  direc¬ 
tions  and  upwind-differencing  in  the  (  direction.  The  flux  vector 
F  has  been  split  into  and  F~  according  to  its  eigenvalues 
and  the  algorithm  is  given  by 


[l  +  /i«|(A+)"  +  hS(d^  -  hRe-H^J-^M^J  -  A|<] 
x[l  +  hS((A-)^  +  hS^B’'  -  Al,]  (AQ"  -  AQ:)  = 

-  A<{5j‘[(F+)’‘  -  -I-  «/[(F-)"  -  Frl 

-h  -  5.)  +  -  i.)  -  Re-H((S^)} 

The  factored  operators  are  solved  by  sweeping  in  the  (  direction 
and  inverting  tridiagonal  matrices  with  5x5  blocks  for  the  other 
two  directions.  Currently,  a  significant  part  of  the  computational 
time  is  taken  to  form  the  plus  and  minus  Jacobian  matrices  for 
the  flux  vector  F  with  this  numerical  scheme.  Artificial  dis¬ 
sipation  terms  (second-  and  fourth-order)  have  been  added  in 
the  central  differencing  directions  for  stability  reasons.**  The 
numerical  code  is  vectorized  for  the  Cray-2  supercomputer. 

A  body  conforming  finite-difference  grid  has  been  used  for 
the  rectangular  blade  having  an  aspect  ratio  of  7  and  consists 
of  warped  spherical  0-0  grid  topology.  The  flowfield  grid  is  nu¬ 
merically  generated  using  the  three-dimensional  hyperbolic  grid 
solver  of  Steger  and  Chaussee**  with  proper  clustering  in  the 
leading  and  trailing  edge  regions  and  in  the  tip  region.  The  grid 
is  nearly  orthogonal  at  the  surface  and  the  spacing  in  the  nor¬ 
mal  direction  at  the  surface  is  chosen  to  be  0.00006C.  Although 
most  of  the  computations  were  done  with  a  grid  topology  having 
155  points  in  the  periodic  direction  around  the  airfoil,  and  66 
points  each  in  the  spanwise  and  normal  directions,  for  a  total  of 
about  700,000  grid  points,  some  calculations  were  also  done  on 
a  smaller  grid  Jl  x  99  x  31  .  The  outer  grid  boundary  is  chosen 
to  be  at  10  chords  in  all  directions  for  both  sizes. 
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The  boundary  conditions  consist  of  surface  boundary  con¬ 
ditions  and  fariield  boundary  conditions  and  are  applied  explic¬ 
itly.  At  the  rotating  blade  surface  because  of  the  noslip  condi¬ 
tion,  the  contravaiiant  velocities  U,  V  wd  W  are  set  to  zero. 
However,  the  time  metrics  (t,  and  are  nonzero  as  the  blade 
(and  the  grid  attached  to  it)  is  moving  in  azimuth.  The  den¬ 
sity  at  the  wall  is  determined  by  assuming  an  adiabatic  wail 
condition.  The  pressure  along  the  body  surface  is  calcidated 
from  the  normal  momentum  relation  (see  for  example  Ref.  11). 
Having  calculated  the  density  and  pressure,  the  total  energy  is 
determined  from  the  equation  of  state. 

At  the  farfield  boundary  the  flow  quantities  are  either  fixed 
or  extrapolated  from  the  interior  depending  on  whether  the  flow 
is  subsonic  or  supersonic  and  if  it  is  of  inflow-  or  outflow-type  at 
the  boundary.  The  chuacteristic  velocities  of  the  Euler  equa¬ 
tions  determine  the  number  of  flow  properties  to  be  specified 
to  control  the  reflections  of  waves  from  the  boundaries.  For 
the  subsonic-inflow  boundary,  four  quantities  must  be  speci¬ 
fied.  Thus,  density  is  extrapolated  while  the  velocities  and 
the  total  energy  are  specified  by  the  free  stream  values.  For 
the  supersonic-inflow,  all  flow  quantities  are  specified.  At  the 
subsonic-outflow  boundaries,  only  one  quantity  is  specified,  viz., 
pressure  is  fixed.  For  the  supersonic-outflow  condition  all  flow 
quantities  are  extrapolated  from  the  interior.  At  the  plane  con¬ 
taining  the  blade  root  the  condition  dQ/dy  =  0  is  imposed. 
The  interacting  vortex  is  assumed  to  have  an  analytical  repre¬ 
sentation  and  is  initialized  as  described  in  Ref.  7.  The  induced 
velocity  due  to  this  line  vortex  is  caJculated  using  Biot-Savart’s 
Law. 

RESULTS 

The  numerical  code  modified  for  the  purpose  of  calculating 
the  blade-vortex  interaction  has  been  validated  agaiinst  experi¬ 
mental  data  of  Catfadonna  and  Tung*’  for  nonlifting  and  lifting 
hover  conditions  at  both  subcritical  and  supercriticad  tip  speeds 
in  Ref.  10.  The  boundairy  layer  is  assumed  to  be  turbulent  in  all 
catses  considered  here  an  algebraic  eddy  viscosity  model'* 
is  used  for  calculating  the  turbulent  viscosity.  Figure  3  shows  a 
typical  result  for  a  lifting  hovering  rotor  from  Ref.  10;  the  com¬ 
parison  of  the  calculations  with  the  experimental  data  is  very 
good. 

Before  considering  the  interaction  of  the  rotor  flowfield 
with  a  line  vortex  convecting  with  the  free  stream,  the  flow- 
field  of  a  rotor  in  forward  flight  has  to  be  computed.  As  noted 
before*,  at  the  subcritical  flow  conditions  corresponding  to  a 
Mup  =  0.6,  the  time  lag  effects  and  the  three-dimensionad  effects 
appear  negligible  as  shown  in  Fig.  4  and  this  is  the  reason  that 
two-  dimensional  calculations  gave  excellent  agreement  with  ex¬ 
periments  of  Caradonna  et  a!.‘  not  only  for  the  basic  rotor  flow 
in  the  absence  of  vortex  interaction,  ?s  seen  in  Fig.  4,  but  also 
for  the  blade-vortex  interaction  flowfield  as  shown  in  Fig.  2.  It 


should  be  noted  that  the  experiments  of  parallel  blade-vortex 
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Fig.  3  Surface  pressures  of  a  lifting  rotor  in  hover.  Mup  =  0.877, 
e  =  8.0,  Re  =  3.96  x  10*. 


Fig.  4  Instantaneous  surface  pressures  for  a  two-bladed  rotor. 

Mtip  =  0.6,  ft  =  0.2. 

interaction  of  Caradonna  et  al.  were  presented  as  test  case  re¬ 
sults  for  validating  the  two-dimensional  numerical  calculation 
procedures.  Although  this  data  would  provide  a  data  source  for 
validating  the  two-dimensional  codes  at  the  subcritical  flow  con¬ 
ditions,  this  may  not  be  true  for  supercritical  flow  conditions  as 
demonstrated  by  Srinivasan  et  al.* 

At  the  supercritical  tip  flow  conditions  corresponding  to 
Miip  =  0.8  with  an  advance  ratio  of  0.2,  the  basic  rotor  flow- 
field  is  dominated  by  the  shock  wave  and  earlier  two-dimensional 
calculations  showed  strong  flow  three-dimensionality  and  un- 
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steady  time-lags  in  shock  formation  and  decay.  Therefore,  two- 
dimensional  calculations  failed  to  calculate  accurately  even  the 
basic  rotor  flow.  As  shown  in  Figs.  5-6,  the  agreement  of  two- 
dimensional  calculations  with  those  of  three-dimensional  calcu¬ 
lations  and  of  experiment  appears  excellent  in  only  the  first 
quadrant  of  the  rotor  azimuth.  In  the  begining  of  the  second 
quadrant,  the  shock  wave  gets  very  strong  as  the  relative  blade 
speed  reaches  Mach  one.  For  this  flow,  the  two-dimensional 
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Fig.  5  Two-dimensional  calculations  of  basic  rotor  flow  com¬ 
pared  with  three-dimensional  experimental  data.  Mup  =  0.8,  n 
~  0.2. 


Fig.  0  Three-dimensional  Euler  calculations  (Ref.  17)  Com¬ 
pared  with  experiments  (Ref.  1).  Mup  =  0.8,  /i  =  0.2. 


assumption  appears  inadequate.  The  strengthening  of  the  shock 
wave  tmd  the  eventual  progressive  demise  have  strong  time- 
lags  built  in  to  it.  The  flow  three-dimensionality  appears  to 
be  accentuated  in  the  presence  of  the  shock  wave.  Because  of 
these  features,  the  two-dimensional  results  greatly  overpredict 
the  strength  of  the  shock  wave  even  for  the  basic  rotor  flow  (in 
the  absence  of  a  vortex  interaction).  Of  course,  with  a  vortex 
interaction  the  flowfield  gets  further  complicated.  These  and 
other  important  results  will  be  presented  in  the  final  version  of 
the  paper. 

CONCLUSIONS 

A  numerical  procedure  is  presented  to  calculate  the  tm- 
steady,  viscous,  three-dimensional  flowfields  of  a  helicopter  ro¬ 
tor  in  forward  flight  and  blade-vortex  interactions  in  subcritical 
and  supercritical  flow  conditions.  Important  flow  features  such 
as  time  lag  effects  in  shock  growth  and  demise,  as  well  as  the 
importance  of  three-dimensional  effects  are  discussed.  The  test 
cases  considered  here  for  computations  correspond  to  the  two 
sets  of  test  data  generated  by  Caradonna  et  al.’,2  on  a  model 
two-bladed  rotor  in  a  wind  tunnel. 
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